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The properties of higher-index saddle points have been invoked in recent theories of the dynamics 
of supercooled liquids. Here we examine in detail a mapping of configurations to saddle points 
using minimization of |VZ?| 2 , which has been used in previous work to support these theories. 
The examples we consider are a two-dimensional model energy surface and binary Lennard- 
Jones liquids and solids. A shortcoming of the mapping is its failure to divide the potential 
energy surface into basins of attraction surrounding saddle points, because there are many 
minima of |Vi5| 2 that do not correspond to stationary points of the potential energy. In fact, 
most liquid configurations are mapped to such points for the system we consider. We therefore 
develop an alternative route to investigate higher-index saddle points and obtain near complete 
distributions of saddles for small Lennard-Jones clusters. The distribution of the number of 
stationary points as a function of the index is found to be Gaussian, and the average energy 
increases linearly with saddle point index in agreement with previous results for bulk systems. 
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I. INTRODUCTION 

The properties of the potential energy surface (PES), 
or landscape, of supercooled liquids and glasses have 
been the focus of much recent theoretical research into 
glasses. The origins of this approach date back to Gold- 
stein, who suggested that the dynamics could be sepa- 
rated into vibrational motion about a-.minimum on the 
PES and transitions between minima.tl This idea led to 
the pioneering inherent structure' approach of Stillingcr 
and coworkers BE In this approach the PES is partitioned 
into basins of attraction surrounding the minima (inher- 
ent structures), where a basin of attraction is defined 
as the set of points for which steepest-descent pathways 
lead to the same minimum. This mapping allows a con- 
ceptual framework to be built in which the role of the 
minima can be separated from the effects of vibrational 
motion. The dynamics can then be viewed in terms of the 
transitions between these basins which occur when the 
system passes along a transition state valley. Therefore, 
from the inherent structure point of view the key points 
on the PES are the minima and transition states, which 
are defined as stationary points with Hessian index one, 
i.e. one negative eigenvalue. Higher-index saddle points 
(with I > 2 negative eigenvalues) need not be consid- 
ered in this description, because, by the Murrcll-Laidlcr 
theorem, if two minima are connected by an index two 
saddle, then there must be a lower-energyr-path between 
them involving only true transition statesfl 

This inherent structure approach has provided impor- 
tant insightsJnto the behaviour, of supercooled liauids 
and glasses, lIlI as well as clustersQ and biomoleculesB For 
example, changes in the dynamics of supercooled liquids 



as the temperature is decreased must correspand to de- 
scent down the PES to lower-energy minimaH Further- 
more, a careful investigation of the density dependence of 
the properties of the sampled minima has suggested how 
changes in the PES can lead to a change in the dynam- 
ics from strong to fragile.113 This approach has also been 
applied to investigate non-equilibrium properties: age- 
ing involves a slow decrease in the energy of the sampled 
minima as the system heads towards equilibrium. Ejc£ 

Much of the above work has focussed on the poten- 
tial energy landscape as sampled under particular con- 
ditions of density and temperature. It is also of interest 
to determine the fundamental characteristics of the land- 
scape, which is, of course, independent of temperature, 
atomic masses and coordinate system. For example, one 
of the most important properties is the distribution of 
minima. Exhaustive enumeration of the minima of small 
systemsl3E-3 seems to confirm the theoretical conjecture 
that—the number of minima increases exponentially with 
sizeBE£l The distribution of minima as a function of the 
potential energy can also be obtained by inverting sim- 
ulation data. This iipersion was first performed fejj, 
medium-sized cluster ,&3 and later for model glassestlll 
revealing that the distribution is Gaussian. This-tech- 
niqueJa,as since been applied to supercooled wateiilj and 
silica. E3 More recently, attention has begun to focus on 
the harder task of characterizing the distwisiitions of tran- 
sition states and the resulting barriers 

Alternatives to the inherent structure approach have 
been proposed. In the instantaneous normal mode the- 
ory, developed by Keyes and coworkers^ the focus is 
on the spectrum of Hessian eigenvalues for instantaneous 
configurations. It is argued that diffusion and the asso- 
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ciated barriers are related to the negative eigenvalues of 
this spectrum, altho ugh this idea has been the subject 
of some debateB-aEa-Ea Many of the negative eigenval- 
ues result from the anharmonicity within a well, how- 
ever once these are removed, simulations have indicated 
that the number of diffusive directions in a supercooled 
liquid tends . towards zero near to the mode-coupling 
temperatureE3o 

Another recent proposal emphasizes the role of higher- 
index stationary points on the PES, and attempts to ex- 
plain the origins of strong and fragile liquids in these 
termsEa One of the underlying ideas is that as the size 
of the system becomes large, most of the configuration 
space volume in a basin of attraction is concentrated near 
to the borders of that basin, and so a randomly chosen 
point is more likely to lie closer to a saddle point than 
to a minimum.Llj Therefore, the proposal is to divide the 
potential energy surface into "basins of attraction" that 
surround stationary points of any index. However, with 
the steepest-descent mapping of the inherent structure 
approach, the basins of attraction only surround minima, 
and convergence to a higher-index saddle point can only 
occur when the starting point lies exactly on the bound- 
ary between two basins of attraction. The volume asso- 
ciated with these boundaries is of measure zero. There- 
fore, a different mapping is required to divide configu- 
ration space in the desired way. Borrowing a trick .that, 
has been previously used to locate transition states j£3E3 
a mapping has been suggested in which steepest-descent 
paths on the function |V£?| 2 , the modulttSpof the gradi- 
ent of the energy squared, are followedHEa Stationary 
points on the PES of any index correspond to minima of 
this new function. 

From the above ma^ning-jit was suggested that be- 
low the mode-couplingoOEH temperature, T c , the sys- 
tem samples minima but that above this temperature 
the average saddle point index increases linearly with 



temperature.! 



Such behaviour has been interpreted 



as a transition from dynamics between basins of minima 
to dynamics between basins of saddles. This approach is 
becoming more widely applied. SESlJO'LJ For example, 
it has been used to analyse the dynamics of ageing. Af- 
ter a quench or crunch (a sudden increase in density) to 
a state that lies below T c for that density, initially the 
system is associated with saddle points whose index de- 
creases logarithmically with time until a crossover itip* 
is reached when the system resides near to minima.E3'c2l 
Here we look in more detail at the IVi^ 2 mapping 
and how well it achieves its aim of dividing configura- 
tion space into neighbourhoods around stationary points 
of any index. We first examine in Section || a model 
two-dimensional energy surface that can easily be visu- 
alized. Then, in Section III we study the properties of 



the mapping for a much-studied binary Lennard-Jones 
system. Given the problems with the iVi^ 2 mapping 
we then follow an alternative approach to studying the 
properties of higher-index saddle points. In section IV 



for small Lennard-Jones clusters and then analyse their 
properties. Finally, we conclude with a discussion of some 
of the issues raised by our results in relation to recent 
work. 



II. MULLER-BROWN SURFACE 

We first examine the effect of the |V-E| 2 mapping for a 
model two-dimensional energy surface that we can easily 
visualize. We use the Miillcr-Brown surface,!^ which has 
the form: 



E(x,y) 



where 



4 

£ 

i=l 



exp[a I (2;-2;°) 2 +6 l (x-^)(y-y°)+c l ( 2 ;-y°) 2 ] 

(1) 



A = (-200, -100, -170, 15) a = (-1, -1, -6.5, 0.7) 
b = (0,0,11,0.6) c= (-10,-10,-6.5,0.7) 

x° = (1,0,-0.5,-1) y° = (0,0.5,1.5,1). (2) 

This surface has been used as a test system for local op- 
timization algorithms, jmrl its properties have been thor- 







we obtain near complete distributions of saddle points 



oughly examined £3 

Contour plots of E and |Vi?| 2 are shown in Figure 
|l| and information concerning the minima of |Vi?| 2 are 
given m Table @Ma The section of the energy surface that 
we consider has five stationary points: three minima and 
the two transition states that connect them. Stationary 
points of the |Vi?| 2 surface satisfy 2Hg = 0, where H 
is the Hessian (second derivative matrix) and g = V-B 
is the gradient vector. Obviously, stationary points of 
the PES have iVi^ 2 — 0, and so correspond to minima 
of |V-B| 2 . However, there are additional minima on the 
|V-E| 2 surface with |Vi?| 2 > 0, where g is an eigenvector 
of H with zero eigenvalue, i.e. there is an inflection point 
in the direction of the gradient. We will refer to these 
two types of Vi?| 2 minima by the labels SP (stationary 
point of E) and NSP (non-stationary point of E) . There 
are four such NSP's on the Muller-Brown surface. The 
possibility of NSP's has been previously noted in Refs. 
p9| , |33| and |34||— however it was claimed that their effect 
was negligible .E3 

A further property of the NSP's is that they must lie on 
a gradient extremal, one definition of which is H g = A g. 
Gradient extremals are curves for which each point is an 
cxtremum of \VE\ = |g| along the corresponding energy 
contour. The gradient extremals for the Miiller-Brown 
surface have been calculated in Ref. [l7j 

The basins of attraction associated with the SP's and 
NSP's on the S7E\ 2 surface are shown in Figure |l|. For 
the region of configuration space that we depict here, the 
majority of this space corresponds to basins of attraction 
associated with NSP's. Although much of this configu- 
ration space has a relatively high energy, the NSP basins 
of attraction do extend into some low energy regions. In 
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FIG. 1: Contour diagrams of (a) E and (b) | S7E\ 2 for the 
Miiller-Brown surface. The red lines divide the surfaces 
into basins of attraction surrounding each minimum. In 
(a) the basin boundaries of |V£?| 2 have been superim- 
posed in blue. Points have been added corresponding to 
in (a) the minima of |Vi?| 2 and in (b) the maxima and 
minima of |Vi?| 2 . In (b) contours occur every 1000 in 
the range to 20000, then every 5000 up to 100000 and 
finally every 50000 beyond that range. 



particular, NSP1 is lower in energy than the transition 
state SP5, and on moving away from the minimum SP1 
along the softest mode, the first new basin of attraction 
that is encountered corresponds to NSP1. Furthermore, 
the basins of attraction associated with NSP2 and NSP3 
extend below the energy of SP5 into regions correspond- 
ing to the walls of the basin of attraction associated with 
SP1 on the original surface. Minimization of |Vi?| 2 for 
points from these regions will lead to a considerable in- 
crease in energy. 

It is clear from Figure |l| that the | VE\ 2 surface is much 
more rugged than the original energy surface. Firstly, 
the surface has more minima. Secondly, the ratio of 
maximum to minimum non-zero Hessian eigenvalues of 
the |V.E| 2 function at SP's has a magnitude that is 
roughly the square of the corresponding ratio for station- 



TABLE I: Minima of \VE\ 2 for the Muller-Brown sur- 
face. Those that are also stationary points of E are la- 
belled SP, and those that are not are labelled NSP. The 
index, /, is the number of negative eigenvalues of the 
Hessian at that point. 





I 


E 


|V£| 2 


X 


V 


SP1 





-146.700 


0.0 


-0.558 


1.442 


SP2 





-108.167 


0.0 


0.623 


0.028 


SP3 





-80.768 


0.0 


-0.050 


0.467 


SP4 


1 


-72.249 


0.0 


0.212 


0.293 


SP5 


1 


-40.665 


0.0 


-0.822 


0.624 


NSP1 





-56.235 


10892.8 


-1.169 


0.741 


NSP2 





19.057 


175.2 


-1.559 


1.543 


NSP3 


1 


21.394 


5018.2 


-0.097 


1.076 


NSP4 





27.070 


6533.2 


-0.995 


-0.053 



ary points of E, because the second derivative of IV-EI 2 
includes a product of the original Hessian matrix. For 
SP1 this ratio is ~100 and so the well surrounding SP1 
on the |V.E| 2 surface is extremely asymmetric and nar- 
row. 

This second feature of the | VE\ 2 surface has important 
practical consequences. In the language of optimization 
theory, such points are said to be ill-conditioned,L2l and so 
even simple minimization of |Vi?| 2 is likely to be rather 
slow. These effects will be further examined when we 



consider a binary L J system in Section [II . 

Although the Miiller-Brown surface is not a physical 
example, it does suggest that the division of configuration 
space into basins of attraction surrounding the minima 
of |Vi?| 2 could be problematic. 



III. BULK BINARY LENNARD-JONES 

Binary Lennard-Jones (BLJ) mixtures have been 
extensively studied in an effort to elucidate the complex 
phenomenology of glasses, as they do not crystallize on 

Most of these investigations have reported changes in 
behaviour around the critical., temperature predicted 
by mode-coupling theoryj3ti3'E3 which is calculated as 
T c w 0.435 for the mixture we consider hejje.E3ti3 For 
example, Sastry, Debenedetti and StillingerEl presented 
evidence that non-exponential relaxation starts below 
about T = 1, while the height of the effective barriers 
to relaxation increases significantly around T — 0.45. 
They view the regions below T = 1 and T = 0.45 
as 'landscape-influenced' and 'landscape-dominated' 
regimes, respectively. They also found that local minima 
obtained by quenching from configurations generated 
at T = 0.5 could escape to different local minima 
much more easily than local minima obtained from 
configurations generated at T = 0.4. 



4 



Using an instantaneous normal modes picture Donati, 
Sciortino and Tartaglia also concluded .that activated 
dynamics becomes important around T c .EZl Schr0der et 
al. found that the liquid dynamics could be separated into 
transitions between minima and vibrational motioru at a 
similar temperature,Ej and correlated motions of atoms 
in groups that grew with decreasing temperature were 
reported by some of the same workers. E3 Two different 
groups have recently reported that the typical Hessian in- 
dex of stationary points sampleeLbv BLJ systems extrap- 
olates to zero, again around T c ,BeJ and we will present a 
more detailed investigation of this behaviour below. Our 
results show that transition states are still accessible be- 
low T c , but support the general consensus that the PES 
sampled by the BLJ system changes in character some- 
where around T c . 

The BLJ mixture that we consider involves a 256- 
atom supercell containing 205 (80%) A atoms, and 51 
(20%) B atoms, with parameters oaa — 1-0, <jab = 
(Jbb = 0.88, e A A = 1-00, e AB = 1.5, and e B B = 0.5.E3 
The units of distance and energy were taken as <j aa and 
tAA ■ The chosen box length gives a fixed number density 
of 1.2 a and a cutoff of 2.5 a aa was used along with 
the minimum image convention and a shifting/truncation 
scheme that ensures continuity ftf-the energy and its first 
derivative, as in previous work.oE 2 ] 

Standard microcanonical molecular dynamics (MD) 
simulations of 10 5 equilibration steps, followed by 10 6 
data collection steps, were first run using the Verlet prop- 
agator at a series of increasing total energies (Table |l). 
The starting point for the first run_was the lowest-energy 
minimum found in previous work,Eil and subsequent runs 
used the final configuration of the previous trajectory 
as the starting point. A time step of 0.003 reduced 
units was employed in each case. Every 1000th config- 
uration from the data collection phase was saved and 
used as a starting point for the following geometry op- 
timizations: (1) minimization using a modified version 
of Nocedal's LBFGS algorithm,EJ (2) a tcjwsifion state 
search using hybrid eigenvector-followingJJOEJ (3) miit. 
imization of |Vi?| 2 using Nocedal's LBFGS algorithm.E 3 ] 
The first two searches on the conventional PES employ 
standard techniques!!] and were followed by between one 
to three full eigenvector-following steps to converge the 
root-mean-square (RMS) force below 10~ 7 reduced units 
and check the Hessian index of the stationary point, de- 
fined as the number of negative Hessian eigenvalues. All 
the searches in (1) and (2) are tightly converged to sta- 
tionary points of the required index with the above tol- 
erance, which is more stringent than the criteria used in 
Ref. m 

As pointed out in Section [n] the ratio of the maximum 
to the minimum eigenvalue of | \7E\ 2 near to an SP makes 
minimization of Vi?| 2 laborious. There is a further 
problem because the second derivatives of the shifted- 
truncated BLJ potentials are discontinuous at the cut- 
off, and hence the derivatives of Vi?| 2 have correspond- 
ing discontinuities. Nevertheless, it is possible to con- 



verge the IVi^l minimizations to an root-mean-square 
force below 10 -5 , at which point the value of |VS| 2 is 
generally converged to at least nine decimal places. Such 
accuracy should be acceptable for the present purposes, 
and was achieved by fixing the neighbour list during min- 
imization for SP's and NSP's close to convergence but 
suffering from discontinuities. Typically, these minimiza- 
tions of | Vi?| 2 involve two orders of magnitude more steps 
than minimizations on the original PES. 

None of the-nccvious work reporting results of |V-E| 2 
minimizationsBE.3 has specified any details of the calcu- 
lations, such as the algorithms employed, the convergence 
criteria or the number of saddles actually found. We are 
therefore unable to provide any detailed comparisons. 

The glass transition temperature for this system under 
our simulation conditions lies between 0.4 and 0.5, as is 
evident from the caloric curve (not illustrated) and jumps 
in various quantities tabulated in Tables || and III. Of 
course, the precise temperature at which the glass tran- 
sition occurs depends upon the rate at which the tem- 
perature is changed. The number of different minima 
and transition states located from the 10 3 different start- 
ing points decreases markedly on glass formation. How- 
ever, a significant number of distinct minima and transi- 
tion states are sampled below the glass transition, as ex- 
pected from previous work that revealed large numbepa 
of low barrier "non-diffusive" pathways for this system.Eil 
The fraction of negative eigenvalues, i = I/(3N — 3), 
(the three zero eigenvalues corresponding to translations 
are excluded) located by minimizing IV-EI 2 jumps by a 
factor of about two on melting, and continues to rise 
approximately linearly at higher temperature. This re- 
sult is in. Jine-.with previous calculations for supercooled 
liquids JPJoo and also with studies based upon instan- 
taneous normal mode analysis .EH However, in contrast to 
the claims of Refs. |33|, |j and ^l], the value of i does 
not vanish, even for our low-temperature glassy configu- 
rations. 

Throughout the temperature range studied here the 
vast majority of minimizations of iVi^l 2 converge to 
NSP's. The largest percentage of SP's occurs for run 
2, but this is probably a fluctuation caused by the non- 
ergodic nature of the low temperature simulations. Sim- 
ilarly, in Ref. |34| NSP's were said to be "frequently sam- 
pled" and these points were excluded from the calculated 



properties. By contrast, in Ref. 33 the number of NSP's 



was said to be "negligible with respect to the number of 
true saddles." However, stimulated by the current re- 
sults, the authors of this paper have now identified the 
source of this discrepancy. A reanalysis of their origi- 
nal results has confirmed that the majority rfthe iVi^ 2 
minima that they found were in fact NSP's £3 

For each run we find the mean energy differences be- 
tween the starting point and the structures obtained af- 
ter searching for minima (min) , transition states (ts) and 
minimizing | VE\ 2 (G2) are in the order A£ m i n > A.Et s > 
AEq2- Hence, in terms of energy, the system is closer 
to points found by minimizing \S7E\ 2 than it is to mm- 
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TABLE II: Mean total energy, E, potential energy, PE, kinetic energy, KE, and kinetic equipartition temperature, 
T, for the MD runs. The ± values represent one standard deviation. #min, #ts and #G2 are the number of distinct 
minima, transition states and stationary point of | Vi?| 2 found for 10 3 searches (excluding permutational isomers), 
as described in the text. %SP and %NSP are the percentage of quenches on the |Vi?| 2 surface that converged to 
stationary points and non-stationary points of E, respectively (out of 10 3 total). 



run 


E 


PE 


KE 


T 


^min 


#ts 


#G2 


%SP 


%NSP 








256-atom supercell 












1 


— 1699.879 ± 0.002 


—1753 ± 2 


53 ±2 


0.138 ±0.005 


54 


274 


740 


1.6 


98.4 


2 


— 1599.268 ± 0.005 


— 1704 ± 4 


105 ±4 


0.276 ±0.010 


46 


509 


909 


9.6 


90.4 


3 


— 1500.070 ± 0.008 


—1656 ± 6 


156 ±6 


0.409 ±0.015 


280 


789 


1000 


1.8 


98.2 


4 


-1399.487 ± 0.011 


-1594 ± 7 


194 ±7 


0.510 ±0.019 


987 


1000 


1000 


3.2 


96.8 


5 


-1301.872 ±0.015 


-1540 ±9 


238 ±9 


0.625 ±0.023 


984 


1000 


1000 


3.0 


97.0 


6 


-1201.004 ± 0.020 


-1485 ± 11 


284 ± 11 


0.745 ± 0.028 


991 


1000 


1000 


4.8 


95.2 


7 


-1098.883 ±0.025 


-1431 ± 12 


332 ± 12 


0.871 ± 0.032 


995 


1000 


1000 


5.9 


94.1 


8 


-999.982 ± 0.031 


-1380 ± 14 


380 ± 14 


0.997 ± 0.036 


994 


1000 


1000 


4.8 


95.2 








320-atom crystal 












1 


-2192.679 ± 0.002 


-2270 ± 2 


78 ±2 


0.162 ±0.002 


1 


39 


1 


100.0 


0.0 


2 


-1992.690 ±0.007 


-2175 ±6 


182 ±6 


0.382 ±0.012 


1 


49 


2 


99.4 


0.6 


3 


-1792.710 ±0.013 


-2082 ± 9 


289 ±9 


0.606 ±0.019 


1 


56 


10 


91.0 


9.0 


4 


-1592.739 ±0.020 


-1989 ± 12 


397 ± 12 


0.831 ±0.026 


1 


64 


55 


66.4 


33.6 


5 


-1392.773 ±0.034 


-1889 ± 16 


496 ± 16 


1.041 ± 0.034 


187 


721 


860 


25.5 


74.5 



ima and transition states, especially at high temperature. 
This result is unsurprising. By the same logic as the 
Murrell-Laidler theoremQ one expects the average poten- 
tial energy of a stationary point to increase with /. 

However, the Euclidean distance between the above 
points is practically the same for all three searches (we 
checked that the centre-of-mass did not change during 
geometry optimization). In terms of distance, the sys- 
tem seems to be equally close to a minimum, a true 
transition state, and an SP or NSP at all temperatures. 
Therefore, it cannot simply be claimed that the Vi?| 2 
mapping takes the system to its nearest stationary point. 
These results contrast with the behaviour for a spin glass, 
namely the p-spin spherical model, for which the closest 
minimum is significantly further away from an equilib- 
rium configuration than is the closest saddle when the 
system is above the glass transition.E3 For this model the 
distance to the closest saddle was also found to be inde- 
pendent of temperature, again differing from our results 
for a structural glass (Table III). 

To check that our results do not depend significantly 
on the minimization algorithm employed we repeated 
some of the IVi^l 2 minimizations using a true steepest- 
descent approach. Both fifth-order RungaJCutta and 
Bulirsch-Stoer algorithms were considered ,tiZl with the 
former method proving to be the more efficient, although 
it still required between lOf^-and I0 3 times more steps 
than LBFGS minimization^ To reduce the computa- 
tional cost we used fifth-order Runga-Kutta integration 
of the steepest-descent equations for order. 10 5 steps and 
then switched to LBFGS minimization.!! 3 ] For 100 reg- 
ularly spaced starting points from trajectories 1 and 8 



the statistics produced by this alternative minimization 
scheme were very similar to the results in Tables || and 



[II 



We have also generated results for the BLJ crystal with 
space group. I^/mmm, which we have recently described 
elsewhereS Here the supercell consists of 320 atoms with 
box lengths of 6.1698 (twice) and 7.0053; the other pa- 
rameters are the same as above, with 256 A and 64 B 
atoms. Configurations were saved from five MD runs of 
10 6 steps each, with 10 5 steps of equilibration, as for the 
smaller system. The solid is superheated, and only es- 
capes from the crystal in the highest energy run on this 
time scale (Table ||). The fraction of NSP located in 
Vi?| 2 minimizations increases systematically from zero 
at the lowest energy, where each minimization finds the 
crystal. The / values associated with these runs are all 
much lower than for the smaller supercell, indicating that 
the linear rise in I above a threshold temperature that 
has been observed for supercooled liquids is not simply a 
universal effect of temperature, but is specific to that re- 
gion of configuration space. However, stationary points 
of index up to four are located in run four, where the 
system still does not escape from a single permutational 
isomer of the crystal. On minimizing the energy from 
all the |V£?| 2 stationary points located in the latter run 
97.8% relax to the crystal, but nine other minima are 
also found. This result simply reflects the fact that min- 
imizing |V-E[ 2 can raise the energy sufficiently for the 
configuration to escape from the crystal, even though the 
system is trapped there on a long time scale. 

In view of the above results, and the dominance of 
NSP's on minimizing | V-E| 2 , we question whether the dy- 



6 



TABLE III: Mean energy differences, A, and displacements, D, between the starting point and the converged geometry 
after searching for minima (min), transition states (ts) and minimizing |Vi?| 2 (G2). isp and znsp are the fractions 
of negative Hessian eigenvalues after minimizing |Vi?| 2 , split into stationary points and non-stationary points of E, 
respectively. The ± values represent one standard deviation. 



run 




n 

U mi 11 




As 


AE G2 




JVTQD V 10^ 


7QT=> X 1 










256-atom supercell 








i 
i 


01 ft _|_ 17 
Z40 zt 1 ( 


zz.U zt l.u 


OHO 4-17 

ZUo ± 1 / 


23.1 ± 2.0 


199 ± 16 


Z4.Z rt 1.0 


c 1 i 1 7 
0.1 1 1. / 


1 O I Q 1 

4.z zt Z. 1 


Z 




0^ 1 n s 

ZO.Z zt U.o 


zoo ± zz 


25.4 ± 1.4 


245 ± 22 


O^ (\ 4- 1 O 
ZO.O ± l.Z 


C/j ion 
0.4 zt Z.U 


C I 1 1 Q 

O.l zt l.o 


Q 

o 


oUo ± 04 


o^ A 1 1 7 

Z0.4 zt 1. / 


OOO 4- AO 

Zyy zt OZ 


25.7 ± 1.6 


286 ± 85 


OR 14-17 
ZD. III./ 


c 710 /) 
O. f zt Z.4 


E Q _L 1 7 
O.o zt 1. / 


1 
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0.0 =b 0.0 
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8.8 =b 0.05 


144 ± 10 
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0.0 =b 0.0 


0.0 =b 0.0 
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8.9 ±0.1 
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9.0 ±0.2 
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15.4 =b 3.7 


432 ± 18 
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1.2 zb 1.3 



namics of a supercooled liquid can usefully be described 
in terms of the basins of attraction of iVi^ 2 . Most of 
these basins do not correspond to stationary points of 
E at all for the present system. It is also interesting to 
note that a more recent study instead used a 'Newton' 
method to locate stationary points of any index. cil This 
approach can potentially avoid the NSP problem, how- 
ever few details were provided and the nature of this new 
mapping should also be carefully examined. 

Furthermore, the configuration does not appear to be 
closer (in terms of distance) to stationary points of any 
particular index as the temperature varies. Although our 
results do not support the claim that / vanishes at a well- 
defined temperature above the glass transition, they do 
confirm that there is a dramatic decrease in the num- 
ber of different local minima sampled around this point. 
Since the prediction that / should vanish at Tmct is a 
mean- field result J£j it is perhaps not surprising that it 
is not precisely obeyed. Similarly, the relaxation time 
scale for structural glasses does not actually diverge at 
?mct because activated processes can still occur below 
this temperature. Of course, it should be remembered 
that all our low temperature results for the 256-atom 
cell correspond to non-equilibrium data, since there is a 
crystalline phase available E3 

IV. LENNARD-JONES CLUSTERS 

From simulations one can obtain the probability distri- 
butions of the system being in the basin of attraction of 
a minimum of energy E at a temperature T. As one can 
evaluate the partition function of a minimum within the 



harmonic approximationy-,or more accurately using an- 
harmonic expressions ,r 6 IF°l the probability distributions 
can be i nverted to obtain distributions for the number of 
minimaJlj'PH 18 ! 

However, even if we could find a way of dividing up 
configuration space into basins around saddle points (i.e. 
without the problem of NSP's associated with the | WE\ 2 
mapping), we could not find the actual distributions of 
the number of saddles from simulation. The missing in- 
gredient is an expression for the partition function asso- 
ciated with the basin around a saddlejld without which 
probability distributions of saddles obtained from simu- 
lation cannot be inverted. 

Therefore, an alternative approach is needed to obtain 
distributions of saddle points. Here we aim to obtain 
(near) complete sets of saddle points for a model finite 
system, namely small Lennard- Jones clusters. This task 
has been previously attempted for minima and— transi- 
tion states up to A=13 by Tsai and JordanjU since 
then larger databases have been obtained for some of 
these clusters. fflo Although there are standard meth- 
ods available to find minima and transition states, finding 
saddle points of a particular index represents a new chal- 
lenge. Eigenvector-following provides an efficient way to 
locate transition states, where we search uphili-in one di- 
rection while minimizing in the tangent space.O We have 
simply extended this approach to find a saddle point of 
index / by searching uphill in / orthogonal directions, 
while minimizing in all other directions. All the uphill di- 
rections were treated in the same way as for a transition 
state search, and were orthogonalized to the search direc- 
tion and gradient in the tangent space minimization, ret. 
quiring only minor modifications of our usual approach.!] 
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TABLE IV: Number of saddle points of each index for LJjv clusters. The numbers in italics are likely to be far from 



complete. 
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To generate the sets of stationary points, we begin by 
obtaining samples of minima and transition states as in 
previous applications BEJO From these sets of stationary 
points we search for index two saddles after randomly 
perturbing the coordinates of the minima and transi- 
tion states. We typically perform twenty such searches 
for each stationary point. We then iteratively repeat 
this procedure for higher-index stationary points, at each 
stage performing searches from all stationary points of 
lower index. This procedure is terminated when no sta- 
tionary points of a particular index are found. 

The sets of stationary points obtained in this man- 
ner are typically incomplete, and the incompleteness is 
larger for stationary points of lower index, for which 
fewer searches have been conducted. To converge the 
sets a reverse procedure was performed. Starting from 
the stationary points of highest index (7 ma x) searches 
are performed for saddle points of index I max — 1 fol- 
lowing a random perturbation. Typically five searches 
from each stationary point are enough to ensure conver- 
gence. Searches are then performed for stationary points 
of index J max — 2 from all those of higher index, and 
so on until the searches for minima are completed. The 
importance of this reverse procedure is evident, for exam- 
ple, from the approximately 50% increase in the number 
of LJ13 transition states located when the searches from 
higher-index saddle points were performed. 

The numbers of stationary points as a function of the 
saddle point index are given in Table IV. We were able 
to find complete distributions of stationary points of any 
index for all clusters up to iV=9. However, above this 
size the search had to be terminated at low index, be- 
cause the numbers of stationary points involved are too 
large for a characterization of the whole distribution to 
be feasible. For example, extrapolation suggests that for 
LJ13 that there are of the order of 10 8 saddle points of 
the most common index. For these larger sizes the re- 
verse procedure was still performed, but starting from 
the highest index for which we performed searches. As a 
result the number of saddle points of the highest index 
searched is likely to be much less than the true total, be- 
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FIG. 2: The number of saddle points as a function of 
(a) the saddle point index and (b) the size. In (a) the 
data points are from Table IV, the curves are Gaussian 
fits to the data, each curve is labelled by the cluster size 
and only sizes for which we have obtained complete dis- 
tributions are represented. In (b) lines are labelled by 
the saddle point index. 



TABLE V: Parameters for the Gaussian fits of the dis- 
tributions n S p(I). Results are only included for sizes that 
have a complete distribution. 



N 


max 
"'sp 


-^mid 


2 mid 


a 


4 


1.49 


1.700 


0.283 


1.699 


5 


5.59 


2.989 


0.332 


1.545 


6 


27.9 


4.063 


0.339 


1.663 


7 


191.4 


4.944 


0.330 


1.738 


8 


1706 


5.858 


0.325 


1.791 


9 


18782 


6.728 


0.320 


1.853 



cause no searches from saddle points of higher index have 
been performed. 

One particularly striking feature of the results is the 
large number of higher-index stationary points for sys- 
tems of such small size, especially relative to the num- 
ber of minima. For example, LJg has approximately 800 
times more index seven saddles than minima. 

To quantitatively probe these distributions we plot in 
Figure |^a the number of saddle points against the inten- 
sive measure of the saddle point index, i. For a cluster, 
i = I/(3N — 6) because there are a further three zero 
eigenvalues corresponding to rotations. As is clear from 
Figure ||a the data fits very well to the Gaussian form: 

)• (3) 



n Bp (J) =7i™ ax exp 



2a 2 



mid ) 



The agreement is less good for the smaller sizes, but this 
result is unsurprising considering the small number of 
stationary points. 

The parameters of the Gaussian fits appear in Table [v|. 
It is particularly noteworthy that the mid-points of the 
distributions in Figure ||a are approximately constant. 
Wd ~ 1/3 which implies I m id ~ N ~ 2. Of course, the 
tail of the Gaussian is cut off at 1—0, but it is also cut 
off beyond 2/ m jd. There are no stationary points for / > 
7 max = 2A - 4. 

Another interesting feature of the Gaussian distribu- 
tion is that the standard deviation of the distributions is 
only a weak function of N, scaling sublinearly with size. 
Therefore, the distributions when considered as a func- 
tion of i (rather than I) , as in Figure ^, become narrower 
as the size increases. 

Eq. (||) allows us to predict the ratio of the number of 
transition states to minima: 



exp 



2/mid — 1 



exp 



/2AT- 5 
V 2a 2 



(4) 



This ratio scales less than exponentially with N, because 
a is a weakly increasing function of N. The above equa- 
tion can be rearranged to obtain an expression for a: 



AT -5/2 



log (n ts /n r , 



(5) 




FIG. 3: The ratio of the number of minima to transition 
states as a function of N. 



The value of a thus obtained will, of course, involve a 
greater error than that obtained from fitting to the com- 
plete saddle point distribution, but it can be applied to 
the larger clusters for which we do not have complete 
distributions. Again we find that a continues to increase 
slowly with N. 

In Figure ^|b we show how the number of stationary 
points with a particular index depends upon N. Most of 
the plots seem to tend to straight lines for larger N. Of 
course, there are likely to be significant deviations when 
n sp (I) < 100, where the values could reflect peculiarities 
of these small sizes. Furthermore, there is insufficient 
data to definitely confirm an exponential scaling with 
size. 

There is a theoretical-expectation for n m j n to scale ex- 
ponentially with size.0iia Here we present the simple ar- 
gument of Ref. ||, which applies to a sufficiently large 
system. If a system of mN atoms can be divided into 
m equivalent sub-systems of N atoms, and the stable ar- 
rangements of the subsystem are effectively independent, 
then 

n min (mA) = n min (A) m (6) 
The solution of this equation is 

Tl rain 

(N) = exp(aN). (7) 

A similar argument can be given for the number of tran- 
sition states. Assuming the rearrangements associated 
with the transition states can be localized to one sub- 
cell, the whole mA-atom system will be at a transition 
state when one of the subsystems is at a transition state 
and the rest are at a minimum. Therefore, 

n ts {mN) = mn min (A) m - 1 n ts (A) (8) 

The solution of this equation is 

n ts (N) = Aexp(aA). (9) 



TABLE VI: Numbers of minima, n m i n , and transition 
states, Jits, for M13 clusters at three values of the range 
parameter p. 



p 


3 


4 


6 


10 


7^ min 


7 


159 


1477 


10 814 




47 


1366 


26 431 


> 746 283 


Tits /^min 


6.7 


8.6 


17.9 


> 69.0 



Hence, n ts /n m i n grows linearly with size. Our databases 
are consistent with this trend (Figure ||). The value of 
a is system dependent, and we illustrate this fact for 
the 13-atom cluster bound by the Morse potential, M13, 
as a function of the range parameter, p. The results 

>rk,0 and were 



in Table VI supersede those in previous worl 
obtained using transition state searches based on starting 
points obtained by considering hard sphere collisions, as 
described elsewhere.EJ 

This approach to finding expressions for the distri- 
bution of saddles becomes more problematic for saddle 
points of higher index. As / increases, the equations are 
increasingly hard to solve, as more combinations of sad- 
dle points of different index have to be considered, and 
the assumption of sub-system independence becomes less 
plausible. Therefore we need a different approach if we 
are to justify the Gaussian distribution for n sp (I). 

For any cluster there will always be a single stationary 
point of index 2N— 4 corresponding to a linear chain. For 
this configuration there are five zero Hessian eigenvalues 
(three translations and only two rotations), N— 1 positive 
eigenvalues corresponding to bond stretches, as well as 
the 2N — 4 negative eigenvalues corresponding to bond- 
angle deformations. 

The chain contains the minimum possible number of 
nearest-neighbour contacts for a bound configuration. To 
produce any more negative eigenvalues would require the 
dissociation of an atom, but then the system would no 
longer correspond to a cluster of N atoms. Therefore, the 
linear chain must correspond to the saddle point with the 
highest index, in agreement with the value of / ma x that 
we found for each distribution. In fact the linear chain 
is rather exceptional because of its five zero eigenvalues. 
Any other configuration with A — 1 bonds must be non- 
linear and so has 2A — 5 negative Hessian eigenvalues, 
six zeros, as well as the N — 1 positive eigenvalues. The 
linear chain's five zero eigenvalues allow the possibility of 
another negative eigenvalue. 

The analysis of the linear configuration suggests that 
the stationary point index corresponds to the number of 
bond-angle degrees of freedom that have negative eigen- 
values. The stationary points in Figure |] illustrate this 
trend. Of course, there are some stationary points where 
there are negative eigenvalues in the bond-stretch degrees 
of freedom, but these are in the minority. 

This upper limit to / also suggests that the Gaussian 




FIG. 4: Some LJ7 stationary points. Each structure 
is labelled by the value of /, the index of the stationary 
point, and each is at the midpoint of the energy distri- 
bution for stationary points of that index. 



distribution may be a result of the number of different 
ways of choosing negative eigenvalues for 2N — 5 bond- 
angle degrees of freedom (assuming six zero eigenvalues 
and N — 1 bond stretches). This assumption gives a 
binomial distribution 

nM = ^ N 5 % V (10) 
which can be well approximated by a Gaussian with 

Inud = N - 5/2, 

„ max _ (2N-5)\ 



"sp 



W- 5/2) 1 
N -5/2 



(11) 



This simple analysis gives properties that are in good 
agreement with our actual distributions. I m id almost ex- 
actly matches the observed value and a is a weakly in- 
creasing function of N. The main error is the result that 
"-sp(O) = n sp (2N — 5) = 1. This error occurs because 
even when we have the maximum or minimum number of 
bond-angle degrees of freedom with negative eigenvalues, 
there are still a number of different structural ways that 
this can be achieved (roughly exp(aA) in fact). There- 
fore, to correct for this error, the above expression for 
n™ ax can be multiplied by an exponential function fit- 
ted to the number of minima. Also, the Gaussian fit 
to the binomial distribution gives the erroneous predic- 
tion that rits/rimin is independent of N. However, the 
binomial distribution itself gives the correct result, while 
the Gaussian approximation begins to break down at the 
tails of the distribution. 

One quantity that has been focussed upon in previous 
studies of saddle points in glasses has been the va, 
of the saddle point index with potential energyHEJ 
The averaging performed to obtain this function can be 
done in two ways. One can either look at (E(I\), the 
average energy of saddles with index / (Figure ^a), or 
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FIG. 5: I/3N versus E/N . In (a) the points are averages 
over J and in (b) averages over the energy. In (b) the 
curves are cut off at I ma x — 1/2 for those sizes that have 
incomplete saddle point distributions. 



at {1(E)}, the average index of stationary points with 
energy E (Figure |^b). For N > 10 as the energy in- 
creases the latter function saturates at the highest index 
for which we have performed searches. To avoid this ef- 
fect of the incomplete distributions, we have cut off the 
function at / max — 1/2. 

Both averages show an approximately linear increase 

similar to that ob- 
I However, the slope 



of the energy with saddle pointi 
served for supercooled liquids O 
is somewhat lower for (1(E)), because the energy at 
which saddle points of a particular index are the most 
common does not necessarily correspond to the energy 
at the maximum of their distribution, i.e. (E(I)) (Figure 
||). For / < I m id the majority of the range for which they 
are most numerous has E < (E(I)) and for I > 7 m id has 
E>(E(I)). 

The slopes of curves such as those in Figure o have 
been interpreted in terms of a single characteristic bar- 
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FIG. 6: The number of LJg saddle points of index / as 
a function of energy. Each curve is labelled by the index 
/, except for the sum of all the distributions. 



rier height for the PES.0EH As the lines in Figure || 
have similar slopes this 'barrier' is only a weak function 
of cluster size. However, there is of course a distribution 
of barrier heights.EirE3 Furthermore, the experimentally- 
observed activation energy to structural relaxation is un- 
likely to correspond to a barrier associated with a single 
rearrangement, but to the overall, barrier associated with 
a sequence of rearrangements.QO It is also easy to show 
that the average barrier between minima and transition 
states, (A), is different from the slope of Figure ||a: 

(A} = (E ta )-J2<AJ2n ta (12) 

i 

? (E ts ) - (E min ), (13) 

where E l min is the energy of minimum i and n[ s is the 
number of transition states connected to that minimum. 
(A) is likely to be larger because the average over minima 
in the second term of equation ( |l2| ) is usually weighted 
towards the lower-energy minima, since they are con- 
nected to more transition states. For example, for LJ13 
(A) = 1.771e, whereas {Ets) - (£ m in)=0.728e. 

Another property of saddle points that has received 
attention is the lowest eigenvalue of the Hessian. For a 
binary Lennard-Jones liquid a linear decrease in the av- 
erage value of the lowest Hessian eigenvalue was observed 
as the energy is increased above the threshold energy al 
which higher-index saddle points were first observed.Efl 
This result is equivalent to a linear decrease with i. The 
behaviour of this property for our cluster saddle points is 
depicted in Figure [?]. The curves seem to have a common 
form, in which the lowest eigenvalue reaches a minimum 
close to I m id- However, these values of i are much larger 
than those probed in Ref. |34|, and the initial parts of the 
curves in Figure [?] seem to show greater linearity as the 
size is increased. 
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FIG. 7: The average of the lowest eigenvalue of the 
Hessian for saddle points of the same index as a function 
of i. The different curves correspond to different sizes as 
labelled. The eigenvalues are in units of e/a 2 . 



V. DISCUSSION 

So far we have examined the properties of a particu- 
lar mapping that attempts to provide a definition of the 
neighbourhood of a saddle point, and looked at the prop- 
erties of higher-index saddle points for systems where 
we can obtain complete distributions. Here, we want to 
think more about the role of higher-index saddle points 
in the dynamics, assuming that one could find a mapping 
that divides all of the PES in neighbourhoods around the 
closest saddle. First we review the relevant aspects of 
the inherent structure approach to the dynamics, i.e. in 
terms of the dynamics of transitions between the basins 
of attraction surrounding minima. This .approach can be 
formulated in term of a master equatiorJZa describing the 
evolution of the occupation probability of a particular 
minimum in terms of thepfja^es of probability flow into 
and out of that minimum f 



dPjjt) 
dt 



(14) 



where Pi is the occupation probability of the basin 
around minimum i and kij is the rate constant for a tran- 
sition from basin j to basin i. This set of equations can 
then be solved to give a complete picture of the inter- 
minimum dynamics. 

The only assumption that this approach relies upon 
is the Markovian nature of the underlying dynamics, i.e. 
the probability of the transition i — ► j must not depend 
on the history of reaching minimum i, so that the k^ 
are constants for a given temperature or energy. This 
assumption will certainly be true when states within a 
basin of attraction equilibrate on a time scale faster than 



transitions to different minima, i.e. Ti n t cr 3> Ti n t ra , where 
these two time scales are for interbasin and intrabasin re- 
laxation. Indeed, it is this separation of time scales that 
makes the inherent structure approach a natural way to 
describe the dynamics. The breakdown of the Marko- 
vian character of the interbasin dynamics places an up- 
per bound on the temperature for which this approach 
is applicable. Previous results for small alkali halide and 
LJ clusters show reasonable agreement between MD-caLes 
and model rate theory for relatively high energies .oEil 

In the above theory for isomerization rates there is no 
requirement for the system to lie close to the minimum 
in configuration space. During the occupancy of a given 
catchment basin the system could be found, on average, 
close to the boundary. Therefore, the increase in / with 
temperature seen for supercooled liquids above a thresh- 
old temperature does not necessarily imply that the in- 
herent structure dynamics approach has broken down. 
Rather the test is whether the interbasin dynamics are 
no longer Markovian. 

One of the main advantages of equation ( |l4| ) is that 
we caa-jdraw on the mature field of unimolecular rate 
theoryta to calculate the kij. For example, the clas- 
sical limits for the microcanonical and canonical Rice- 
Ramsperger-Kassel-Marcus rate constants in the har- 
monic approximation areH 



k{E) 



E-E^ 
E 



s-l 



-,8-1 ■ 



k(T) = exp(-E^ /kT)- 



(15) 



where E and E< are the total energy and the potential 
energy of the transition state relative to the energy of 
the minimum, s is the number of vibrational degrees of 
freedom, and v and are the geometric mean vibra- 
tional frequencies of the minimum and transition state, 
respectively. 

Angelani et al. have argued that "diffusion is entropy 
driven, even below Tmct" based on the observation that 
the instantaneous potential energy of the systemJies well 
above the SP's obtained by minimizing |V£| 2 .EI How- 
ever, increasingly large values of the total energy are 
needed to obtain significant rates as the By,mber of de- 
grees of freedom increases (see for example^) , as is clear 
from equations ( |l5| ) . To maintain a finite microcanonical 
rate constant as s rises the ratio E/E^ must increase. 
From the canonical viewpoint, the higher total energy is 
simply required to maintain a constant temperature as 
the system size grows. 

In the microcanonical ensemble the entropy is the ap- 
propriate thermodynamic potential, while in the canon- 
ical ensemble we must instead consider free energy bar- 
riers, although it should be remembered that the two 
ensembles are equivalent in the bulk limit. The relative 
energies of the transition state and minimum, and the 
widths of the minimum and transition state valley con- 
tained in the ratio of frequencies, contribute to the rate 
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constant for an elementary unimolecular isomerization in 
both ensembles. Both are included in the standard rate 
expressions above. However, the relative importance of 
different contributions to k(E) cannot be assessed with- 
out considering how the different terms scale with sys- 
tem size. For the canonical rate constant k(T), where it 
is meaningful to consider separate entropic and energetic 
contributions to the free energy barrier, the total energy 
does not enter. The energy difference, E\ for a particu- 
lar class of rearrangement should be an intensive rather 
than an extensive quantity. 

For normal liquids above the melting point Dzugutov 
has demonstrated a strong empirical correlation between 
the entropy and thjc_di£fusion constant .E3 Other recent 
simulation studiesEiJ — have been interpreted using the 
Adam-Gibbs model £3 where the entropy, or part of it, 
enters in a rather different way through heuristic argu- 
ments. On the other hand, atomic diffusion in solids is 
routinely treated by Vineyard's approachJ£3 which is sim- 
ply an application of conventional unimolecular rate the- 
ory to bulk. There is clearly a pressing need to determine 
the limit of applicability of standard rate theory to su- 
percooled liquids, and to develop alternative approaches 
with a firm microscopic basis where necessary. 

Cavagna has suggested that when r m t e r ~ Ti n t ra it 
is more appropriate to think about the dynamics in 
terms of tmiisitions between the neighbourhoods of sad- 
dle points O As equation (|i4j) can be applied to any par- 
tition of the PES into 'basins', not just those surround- 
ing minima, a similar formalism can be developed as 
long as the inter-saddle dynamics are Markovian. How- 
ever, therein lies the problem. If the residence times in 
the basins of attraction surrounding a minimum are too 
short for equilibrium between the vibrational modes to 
be established, it seems even less likely that the necessary 
separation of time scales will hold for basins surround- 
ing saddle points. By definition, there are forces acting 
to take the system out of such regions. Furthermore, in 
contrast to the case of isomerizations between minima, 
there is no established theory for transition rates between 
saddles. In fact, it is unclear how to calculate the parti- 
tion function for the catchment basin of a saddle, which 
would surely be necessary to evaluate rate constants. It 
is therefore hard to see how this view of the dynamics 
can be put on a quantitative footing. In his contribu- 
tion Cavagna simply speculated that the-kitersaddle rate 
constants would increase as I increasedHj 

One other possible criticism of the 'saddles-ruled' ap- 
proach of Cavagna is that it seems to ignore many of 
the effects of the topography of the PES, and so conflicts 
with much of the recent work-that has emphasised the im- 
portance of this topography.BEIEj For example, Cavagna, 
suggests that the origin of strong and fragile behaviourEll 
simply lies in the value of the thermal energy at the tem- 
perature for which saddles begin to be frequently sampled 
relative to the "barrier" obtained from the slope of the 
line giving the averaged-dependence of the index upon 
the saddle point energyEa 



It is also too simplistic to suggest that the index pi. 
a saddle indicates the 'number of diffusive directions '.E3 
In fact, rearrangements may be both diffusive and non- 
diffusive in nature^! and the character can only be di- 
agnosed by calculating steepest-descent paths connecting 
the relevant stationary points, not from purely local in- 
formation. The non-diffusive transitions typically involve 
an atom moving slightly within the cage of its neighbours. 



VI. CONCLUSIONS 

We have examined in detail the behaviour of the | VE\ 2 
mapping of configurations to saddle points, and find 
that it has a number of shortcomings. We agree with 
Cavagnac3 that the mapping cannot partition the whole 
of the PES into basins surrounding the saddle points, 
as claimed by Angelani et al.x3 In fact, for the systems 
we have studied the vast majority of configuration space 
sampled by the supercooled liquid is mapped on to points 
that are not stationary points of the PES. Furthermore, 
the saddle points obtained by this mapping are no closer 
to the initial configuration than are the nearest transi- 
tion state or minimum. Therefore, the mapping does not 
seem to satisfy the requirement of dividing the PES into 
neighbourhoods around the closest saddle. We have also 
pointed out problems with a number of interpretations 
suggested in previous work. 

A more straightforward way to characterize the prop- 
erties of higher-index saddles is for systems where com- 
plete distributions of saddle points can be obtained. Our 
results for LJ clusters reveal that the distributions are 
a Gaussiaa function of the index, as well as of the 
energy.LDES Gaussian index-distributions have also been 
found for random matrices £3 We have suggested an ex- 
planation for this distribution in terms of the number of 
possible ways of assigning negative eigenvalues to a set 
of bond-angle degrees of freedom. We find an approxi- 
mately linear relationship between the potential energy 
of the saddle point and its index, similar to the results 
for supercooled liquids. 

We do not think that a description of the dynamics in 
terms of inter-saddle, rather than inter-minimum, transi- 
tions will offer any advantages, even if a proper partition- 
ing of the PES into catchment basins of saddles can be 
devised. Although the temperature at which the dynam- 
ics become non-Markovian, and hence the linear master 
equation formalism breaks down, has not yet been es- 
tablished it may well be at temperatures where the time 
scale for transitions between minima is comparable to 
the time for equilibrium to be established between the 
vibrational modes. Below this regime we should be able 
to apply standard unimolecular rate theory, but above it 
any approach focussing on uncorrelated transitions be- 
tween local regions of configuration space is likely Jxtfail 
Other approaches, such as mode-coupling theory J?*' 1 
will then be more appropriate. 
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